#include "blockreader.h"
#define VECPRODUCT(a, b) Vec(a.x*b.x, a.y*b.y, a.z*b.z)

int
powerOf2(int val)
{
  int p, q;
  p = 1;
  q = 0;
  while (val > p)
    {
      p*=2;
      q++;
    }
  return q;
}

int BlockReader::minLod() { return m_minLevel; }

int* BlockReader::get1dHistogram() { return m_1dHistogram[0]; }
int* BlockReader::get2dHistogram() { return m_2dHistogram[0]; }

int* BlockReader::get1dHistogram(int nf) { return m_1dHistogram[nf]; }
int* BlockReader::get2dHistogram(int nf) { return m_2dHistogram[nf]; }

uchar* BlockReader::blockLevel() { return m_blockLevel; }

void BlockReader::setTextureSize(int tsz) { m_maxtexSize = tsz; }

void BlockReader::setUserLimits(Vec bmin, Vec bmax)
{
  m_boxMinUser = Vec(bmin.z, bmin.y, bmin.x);
  m_boxMaxUser = Vec(bmax.z, bmax.y, bmax.x);
}

void BlockReader::setMemorySize(int m)
{
  m_memorySize = m;
  int bsz = m_blockSize+m_overlap;
  float frc = (1024.0f*1024.0f)/(m_nVolumes*bsz*bsz*bsz);
  m_maxCacheSize = m_memorySize*frc; // blocks
}

void
BlockReader::init()
{
  m_nVolumes = 1;
  m_currLod = -1;

  for(int nf=0; nf<4; nf++)
    m_hdf5file[nf] = 0;

  m_fileName.clear();

  m_maxCacheSize = 1000; // blocks
  m_blockCache.clear();
  m_lru.clear();

  m_tmpV = 0;

  m_data = 0;
  m_auxdata = 0;
  m_blockLevel = 0;
  m_savedblockLevel = 0;
  for(int nf=0; nf<4; nf++)
    for(int i=0; i<10; i++)
      m_limits[nf][i] = 0;

  m_blockSize = 28;
  m_overlap = 4;
  m_odepth = m_orows = m_ocols = 0;
  m_depth = m_rows = m_cols = 0;
  m_dblocks = m_rblocks = m_cblocks = 0;

  m_maxtexSize = 512;

  m_minLevel = 4;

  setMemorySize(500);

  m_isoval = -1;

  m_summedAreaTable = 0;

  for(int nf=0; nf<4; nf++)
    m_1dHistogram[nf] = 0;
  for(int nf=0; nf<4; nf++)
    m_2dHistogram[nf] = 0;
}

BlockReader::~BlockReader()
{
  for(int nf=0; nf<4; nf++)
    if (m_hdf5file[nf])
      m_hdf5file[nf]->close();
  
  m_fileName.clear();

  if (m_data) delete [] m_data;
  if (m_auxdata) delete [] m_auxdata;
  if (m_blockLevel) delete [] m_blockLevel;
  if (m_savedblockLevel) delete [] m_savedblockLevel;

  for(int nf=0; nf<4; nf++)
    for(int i=0; i<10; i++)
      if (m_limits[nf][i]) delete [] m_limits[nf][i];
  for(int nf=0; nf<4; nf++)
    for(int i=0; i<10; i++)
      m_limits[nf][i] = 0;
  
  m_data = 0;
  m_auxdata = 0;
  m_blockLevel = 0;
  m_savedblockLevel = 0;

  m_odepth = m_orows = m_ocols = 0;
  m_depth = m_rows = m_cols = 0;
  m_compactdepth = m_compactrows = m_compactcols = 0;
  m_dblocks = m_rblocks = m_cblocks = 0;

  QList<qint64> bno = m_blockCache.keys();
  for(int i=0; i<bno.count(); i++)
    delete [] m_blockCache[bno[i]];

  m_blockCache.clear();
  m_lru.clear();

  for(int nf=0; nf<4; nf++)
    if (m_1dHistogram[nf]) delete [] m_1dHistogram[nf];
  for(int nf=0; nf<4; nf++)
    if (m_2dHistogram[nf]) delete [] m_2dHistogram[nf];
  for(int nf=0; nf<4; nf++)
    m_1dHistogram[nf] = 0;
  for(int nf=0; nf<4; nf++)
    m_2dHistogram[nf] = 0;

  if (m_summedAreaTable)
    delete [] m_summedAreaTable;
  m_summedAreaTable = 0;
}

void BlockReader::setIsovalue(int v) { m_isoval = v; }
uchar* BlockReader::blockedData() { return m_data; }
uchar* BlockReader::auxData() { return m_auxdata; }
int BlockReader::blockSize() { return m_blockSize; }

void
BlockReader::gridSize(int &d, int &r, int &c)
{
  d = m_depth;
  r = m_rows;
  c = m_cols;
//  d = m_odepth;
//  r = m_orows;
//  c = m_ocols;
}

void
BlockReader::gridSize(int l, int &d, int &r, int &c)
{
  QString attname = QString("/lod-%1/gridsize").arg(l);
  Attribute attrib = m_hdf5Dataset[0][l].openAttribute(attname.toAscii().data());
  int griddim[3];
  attrib.read(PredType::NATIVE_INT, griddim);
  d = griddim[0];
  r = griddim[1];
  c = griddim[2];
}

void
BlockReader::p2gridSize(int &d, int &r, int &c)
{
  d = m_p2depth;
  r = m_p2rows;
  c = m_p2cols;
}

void
BlockReader::compactGridSize(int &d, int &r, int &c)
{
  d = m_compactdepth;
  r = m_compactrows;
  c = m_compactcols;
}

void
BlockReader::blockedGridSize(int &d, int &r, int &c)
{
  d = m_dblocks;
  r = m_rblocks;
  c = m_cblocks;
}

void
BlockReader::p2blockedGridSize(int &d, int &r, int &c)
{
  d = m_p2dblocks;
  r = m_p2rblocks;
  c = m_p2cblocks;
}

void
BlockReader::setGridSize(int d, int r, int c)
{
  m_depth= d;
  m_rows = r;
  m_cols = c;

  m_dblocks = m_depth/m_blockSize+ (m_depth%m_blockSize > 0);
  m_rblocks = m_rows/m_blockSize + (m_rows%m_blockSize > 0);
  m_cblocks = m_cols/m_blockSize + (m_cols%m_blockSize > 0);

  m_p2depth= qPow(2, powerOf2(m_depth));
  m_p2rows = qPow(2, powerOf2(m_rows));
  m_p2cols = qPow(2, powerOf2(m_cols));

  m_p2dblocks = m_p2depth/m_blockSize+(m_p2depth%m_blockSize > 0);
  m_p2rblocks = m_p2rows/m_blockSize + (m_p2rows%m_blockSize > 0);
  m_p2cblocks = m_p2cols/m_blockSize + (m_p2cols%m_blockSize > 0);
  m_p2dblocks = qPow(2, powerOf2(m_p2dblocks));
  m_p2rblocks = qPow(2, powerOf2(m_p2rblocks));
  m_p2cblocks = qPow(2, powerOf2(m_p2cblocks));


  m_nblocks = m_dblocks*m_rblocks*m_cblocks;
}

void
BlockReader::setFilename(QList<QString> flnm)
{
  m_nVolumes = flnm.count();
  m_fileName = flnm;

  for(int nf=0; nf<m_nVolumes; nf++)
    m_hdf5file[nf] = new H5File(m_fileName[nf].toAscii().data(),
				H5F_ACC_RDONLY);


  for(int nf=0; nf<m_nVolumes; nf++)
    {
      for(int l=0; l<=m_minLevel; l++)
	{
	  QString groupname = QString("/lod-%1").arg(l);
	  try
	    {
	      m_group[nf][l] = new Group(m_hdf5file[nf]->openGroup(groupname.toAscii().data()));
	      QString dataname;
	      dataname = groupname+"/data";
	      m_hdf5Dataset[nf][l] = m_group[nf][l]->openDataSet(dataname.toAscii().data());
	      dataname = groupname+"/minmax";
	      m_minmaxDataset[nf][l] = m_group[nf][l]->openDataSet(dataname.toAscii().data());
	    }
	  catch (FileIException error)
	    {
	      m_minLevel = l-1;
	      break;
	    }
	}
    }

  m_dataType.copy(m_hdf5Dataset[0][0]);
  
  //--------------------------------------
  // allocate memory so that finest level can be loaded
  int d, r, c;
  gridSize(0, d, r, c);

  m_boxMinUser = Vec(0,0,0);
  m_boxMaxUser = Vec(c, r, d) - Vec(1,1,1);

  hsize_t dims[10];
  DSetCreatPropList propList = m_hdf5Dataset[0][0].getCreatePlist();
  propList.getChunk(3, dims);
  int bsz = dims[0]-m_overlap;
      
  int bd = d/bsz + (d%bsz > 0);
  int br = r/bsz + (r%bsz > 0);
  int bc = c/bsz + (c%bsz > 0);
  m_nblocks0 = bd*br*bc;

  if (m_blockLevel) delete [] m_blockLevel;
  m_blockLevel = new uchar[m_nblocks0];
  memset(m_blockLevel, 255, m_nblocks0);

  if (m_savedblockLevel) delete [] m_savedblockLevel;
  m_savedblockLevel = new uchar[m_nblocks0];
  memset(m_savedblockLevel, 255, m_nblocks0);
	  
  if (m_auxdata) delete [] m_auxdata;
  m_auxdata = new uchar[4*m_nblocks0];

  //--------------------------------------
  for(int nf=0; nf<m_nVolumes; nf++)
    {
      for(int lod=0; lod<10; lod++)
	if (m_limits[nf][lod]) delete [] m_limits[nf][lod];
    }

  for(int nf=0; nf<m_nVolumes; nf++)
    {
      for(int lod=0; lod<=m_minLevel; lod++)
	{
	  int d, r, c;
	  gridSize(lod, d, r, c);
	  
	  int bd = d/bsz + (d%bsz > 0);
	  int br = r/bsz + (r%bsz > 0);
	  int bc = c/bsz + (c%bsz > 0);
	  int nblocks = bd*br*bc;
	  m_limits[nf][lod] = new uchar[2*nblocks];
	  
	  hsize_t dims;
	  dims = 2*nblocks;
	  DataSpace ldataspace(1, &dims);
	  DataSpace lmemspace(1, &dims);
	  m_minmaxDataset[nf][lod].read(m_limits[nf][lod],
					m_dataType,
					lmemspace,
					ldataspace );
	}
    }
  //--------------------------------------
      
  //--------------------------------------
  // histogram generation using (possibly) smaller volume
  int ld = 0;
  while((m_nVolumes*m_nblocks0)/qPow(2, ld) > 10000)
    ld++;
  ld = qMin(m_minLevel, ld);

  for(int nf=0; nf<m_nVolumes; nf++)
    generateHistograms(nf, ld);
  //--------------------------------------

  m_summedAreaTable = new ushort[257*m_nVolumes];
}

void
BlockReader::setLodLevel(int slevel)
{
  if (m_currLod == slevel)
    return;

  m_currLod = slevel;

  hsize_t dims[10];

  DSetCreatPropList propList = m_hdf5Dataset[0][m_currLod].getCreatePlist();
  propList.getChunk(3, dims);
  m_blockSize = dims[0]-m_overlap;

  DataSpace dataspace = m_hdf5Dataset[0][m_currLod].getSpace();
  dataspace.getSimpleExtentDims(dims, NULL);
  m_odepth= dims[0];
  m_orows = dims[1];
  m_ocols = dims[2];

  QString attname = QString("/lod-%1/gridsize").arg(m_currLod);
  Attribute attrib = m_hdf5Dataset[0][m_currLod].openAttribute(attname.toAscii().data());
  int griddim[3];
  attrib.read(PredType::NATIVE_INT, griddim);
  setGridSize(griddim[0], griddim[1], griddim[2]);

  //------------------------------
  int cbsz = m_blockSize+m_overlap;
  if (m_tmpV) delete [] m_tmpV;
  m_tmpV = new uchar[m_nVolumes*cbsz*cbsz*cbsz];
}


int
BlockReader::minLod(int isoval)
{
  uchar *vbl = new uchar[m_nblocks0];

  int bsz = m_blockSize+m_overlap;
  for(int lod=0; lod<m_minLevel; lod++)
    {
      int d, r, c;
      gridSize(lod, d, r, c); 
      int bd = d/bsz + (d%bsz > 0);
      int br = r/bsz + (r%bsz > 0);
      int bc = c/bsz + (c%bsz > 0);
      int nblocks = bd*br*bc;
      memset(vbl, 255, nblocks);
      for(int nf=0; nf<m_nVolumes; nf++)
	{
	  uchar *limits = m_limits[nf][lod];
	  for(int bidx=0; bidx<nblocks; bidx++)
	    {	
	      int minval = limits[2*bidx];
	      int maxval = limits[2*bidx+1];
	      if (minval <= isoval && maxval >= isoval)
		{
		  vbl[bidx] = 0;
		  
		  // take neighbouring blocks
		  int ib,jb,kb;
		  ib = bidx/(m_rblocks*m_cblocks);
		  jb = (bidx - ib*m_rblocks*m_cblocks)/m_cblocks;
		  kb = bidx - ib*m_rblocks*m_cblocks - jb*m_cblocks;
		  for(int i=qMax(0,ib-1); i<=qMin(m_dblocks-1, ib+1); i++)
		    for(int j=qMax(0,jb-1); j<=qMin(m_rblocks-1, jb+1); j++)
		      for(int k=qMax(0,kb-1); k<=qMin(m_cblocks-1, kb+1); k++)
			{
			  int idx = i*m_rblocks*m_cblocks + j*m_cblocks + k;
			  vbl[idx] = 0;
			}
		}
	    }
	}

      //---- count valid blocks ----
      int nb = 0;
      for(int bidx=0; bidx<nblocks; bidx++)
	if (vbl[bidx] == 0)
	  nb++;

      int cbsz = m_blockSize+m_overlap;
      int maxb = m_maxtexSize/cbsz;
      if (nb < maxb*maxb*maxb)
	{
	  delete [] vbl;
	  return lod;
	}
    }
  
  delete [] vbl;
  return m_minLevel;
}

int
BlockReader::minLod(uchar *lut)
{
  uchar *vbl = new uchar[m_nblocks0];

  memset(m_summedAreaTable, 0, sizeof(ushort)*m_nVolumes*257);
  for(int nf=0; nf<m_nVolumes; nf++)
    {
      m_summedAreaTable[nf*257 + 0] = 0;
      for(int im=0; im<256; im++)
	{
	  int idx = nf*256*256*4 + 4*im+3;
	  m_summedAreaTable[nf*257 + im+1] = m_summedAreaTable[nf*257 + im] + (lut[idx]>0);
	}
    }

  int bsz = m_blockSize+m_overlap;
  for(int lod=0; lod<m_minLevel; lod++)
    {
      int d, r, c;
      gridSize(lod, d, r, c); 
      int bd = d/bsz + (d%bsz > 0);
      int br = r/bsz + (r%bsz > 0);
      int bc = c/bsz + (c%bsz > 0);
      int nblocks = bd*br*bc;
      memset(vbl, 255, nblocks);
      for(int nf=0; nf<m_nVolumes; nf++)
	{
	  uchar *limits = m_limits[nf][lod];
	  for(int bidx=0; bidx<nblocks; bidx++)
	    {	
	      int minval = limits[2*bidx];
	      int maxval = limits[2*bidx+1];
	      if (m_summedAreaTable[nf*257 + minval] < m_summedAreaTable[nf*257 + maxval+1])
		{
		  vbl[bidx] = 0;
		  
		  // take neighbouring blocks
		  int ib,jb,kb;
		  ib = bidx/(m_rblocks*m_cblocks);
		  jb = (bidx - ib*m_rblocks*m_cblocks)/m_cblocks;
		  kb = bidx - ib*m_rblocks*m_cblocks - jb*m_cblocks;
		  for(int i=qMax(0,ib-1); i<=qMin(m_dblocks-1, ib+1); i++)
		    for(int j=qMax(0,jb-1); j<=qMin(m_rblocks-1, jb+1); j++)
		      for(int k=qMax(0,kb-1); k<=qMin(m_cblocks-1, kb+1); k++)
			{
			  int idx = i*m_rblocks*m_cblocks + j*m_cblocks + k;
			  vbl[idx] = 0;
			}
		}
	    }
	}

      //---- count valid blocks ----
      int nb = 0;
      for(int bidx=0; bidx<nblocks; bidx++)
	if (vbl[bidx] == 0)
	  nb++;

      int cbsz = m_blockSize+m_overlap;
      int maxb = m_maxtexSize/cbsz;
      if (nb < maxb*maxb*maxb)
	{
	  delete [] vbl;
	  return lod;
	}
    }

  delete [] vbl;
  return m_minLevel;
}

bool
BlockReader::initIdentifyValidBlocks(int isoval)
{
  memset(m_savedblockLevel, 255, m_nblocks);

  int d, r, c;
  gridSize(0, d, r, c);
  Vec scale = Vec((float)c/(float)m_cols,
		  (float)r/(float)m_rows,
		  (float)d/(float)m_depth);
  

  //---- set valid blocks ----
  for(int nf=0; nf<m_nVolumes; nf++)
    {
      uchar *limits = m_limits[nf][m_currLod];
      for(int bidx=0; bidx<m_nblocks; bidx++)
	{	
	  int ib,jb,kb;
	  ib = bidx/(m_rblocks*m_cblocks);
	  jb = (bidx - ib*m_rblocks*m_cblocks)/m_cblocks;
	  kb = bidx - ib*m_rblocks*m_cblocks - jb*m_cblocks;
	  
	  Vec bs = Vec(ib, jb, kb)*m_blockSize;
	  Vec be = Vec(ib+1, jb+1, kb+1)*m_blockSize;
	  bs = VECPRODUCT(bs, scale);
	  be = VECPRODUCT(be, scale);
	  bool ok = true;
	  for(int t=0; t<3; t++)
	    ok &= !((be[t] <= m_boxMinUser[t]) || bs[t] >= m_boxMaxUser[t]);
	  
	  int minval = limits[2*bidx];
	  int maxval = limits[2*bidx+1];
	  if (ok && minval <= isoval && maxval >= isoval)
	    {
	      m_savedblockLevel[bidx] = 0;
	      
	      // take neighbouring blocks
	      for(int i=qMax(0,ib-1); i<=qMin(m_dblocks-1, ib+1); i++)
		for(int j=qMax(0,jb-1); j<=qMin(m_rblocks-1, jb+1); j++)
		  for(int k=qMax(0,kb-1); k<=qMin(m_cblocks-1, kb+1); k++)
		    {
		      int idx = i*m_rblocks*m_cblocks + j*m_cblocks + k;
		      m_savedblockLevel[idx] = 0;
		    }
	    }
	}
    }

  //---- count valid blocks ----
  int nb = 0;
  for(int bidx=0; bidx<m_nblocks; bidx++)
    if (m_savedblockLevel[bidx] == 0)
      nb++;
  //----------------------------

  memcpy(m_blockLevel, m_savedblockLevel, m_nblocks);

  int cbsz = m_blockSize+m_overlap;
  int maxb = m_maxtexSize/cbsz;
  bool allblocks = false;
  if (nb < maxb*maxb*maxb)
    allblocks = true;

  return allblocks;

  //QMessageBox::information(0, "", QString("%1 : %2").arg(isoval).arg(nb));
}

bool
BlockReader::initIdentifyValidBlocks(uchar *lut)
{
  int d, r, c;
  gridSize(0, d, r, c);
  Vec scale = Vec((float)c/(float)m_cols,
		  (float)r/(float)m_rows,
		  (float)d/(float)m_depth);
  
  memset(m_summedAreaTable, 0, sizeof(ushort)*m_nVolumes*257);
  for(int nf=0; nf<m_nVolumes; nf++)
    {
      m_summedAreaTable[nf*257 + 0] = 0;
      for(int im=0; im<256; im++)
	{
	  int idx = nf*256*256*4 + 4*im+3;
	  m_summedAreaTable[nf*257 + im+1] = m_summedAreaTable[nf*257 + im] + (lut[idx]>0);
	}
    }

  memset(m_savedblockLevel, 255, m_nblocks);
  //---- set valid blocks ----
  for(int nf=0; nf<m_nVolumes; nf++)
    {
      uchar *limits = m_limits[nf][m_currLod];
      for(int bidx=0; bidx<m_nblocks; bidx++)
	{	
	  int ib,jb,kb;
	  ib = bidx/(m_rblocks*m_cblocks);
	  jb = (bidx - ib*m_rblocks*m_cblocks)/m_cblocks;
	  kb = bidx - ib*m_rblocks*m_cblocks - jb*m_cblocks;

	  Vec bs = Vec(ib, jb, kb)*m_blockSize;
	  Vec be = Vec(ib+1, jb+1, kb+1)*m_blockSize;
	  bs = VECPRODUCT(bs, scale);
	  be = VECPRODUCT(be, scale);
	  bool ok = true;
	  for(int t=0; t<3; t++)
	    ok &= !((be[t] <= m_boxMinUser[t]) || bs[t] >= m_boxMaxUser[t]);

	  int minval = limits[2*bidx];
	  int maxval = limits[2*bidx+1];
	  if (ok && m_summedAreaTable[nf*257 + minval] < m_summedAreaTable[nf*257 + maxval+1])
	    {
	      m_savedblockLevel[bidx] = 0;
	      
	      // take neighbouring blocks
	      for(int i=qMax(0,ib-1); i<=qMin(m_dblocks-1, ib+1); i++)
		for(int j=qMax(0,jb-1); j<=qMin(m_rblocks-1, jb+1); j++)
		  for(int k=qMax(0,kb-1); k<=qMin(m_cblocks-1, kb+1); k++)
		    {
		      int idx = i*m_rblocks*m_cblocks + j*m_cblocks + k;
		      m_savedblockLevel[idx] = 0;
		    }
	    }
	}
    }
  
  //---- count valid blocks ----
  int nb = 0;
  for(int bidx=0; bidx<m_nblocks; bidx++)
    if (m_savedblockLevel[bidx] == 0)
      nb++;
  //----------------------------

  memcpy(m_blockLevel, m_savedblockLevel, m_nblocks);


  int cbsz = m_blockSize+m_overlap;
  int maxb = m_maxtexSize/cbsz;
  bool allblocks = false;
  if (nb < maxb*maxb*maxb)
    allblocks = true;

  return allblocks;
}

void
BlockReader::loadBlocks(QList<Vec> blocklist,
			int slevel)
{
  if (m_currLod != slevel)
    {
      setLodLevel(slevel);
      int bsz = m_blockSize+m_overlap;
      float frc = (1024.0f*1024.0f)/(m_nVolumes*bsz*bsz*bsz);
      m_maxCacheSize = m_memorySize*frc; // blocks
    }
  memset(m_auxdata, 0, 4*m_nblocks);
  memcpy(m_blockLevel, m_savedblockLevel, m_nblocks);

  int nb = 0;
  for(int b=0; b<m_nblocks; b++)
    if (m_savedblockLevel[b] == 0)
      {
	nb++;
	m_blockLevel[b] = slevel;
	m_auxdata[4*b+3] = 100; // to signify that the block is valid
      }

  m_maxlodLevel = slevel;

  int cbsz = m_blockSize+m_overlap;
  int maxb = m_maxtexSize/cbsz;
  int bcount = blocklist.count();
  int cblk = qMin(maxb, bcount);
  int rblk = qMin(maxb, qMax(1, bcount/cblk + (bcount%cblk>0)));
  int dblk = qMin(maxb, qMax(1, bcount/(cblk*rblk) + (bcount%(cblk*rblk) > 0)));

  //QMessageBox::information(0, QString("bcount : %1").arg(bcount), QString("bcount : %1").arg(bcount));

  if (bcount > dblk*rblk*cblk)
    {
      QString flnm = m_fileName[0] + ".error";
      
      {
	QFile file(flnm);
	if (!file.open(QIODevice::WriteOnly |
		     QIODevice::Text))
	  {
	    QMessageBox::critical(0, "Error", "Cannot write to file");
	    return;
	  }

	QTextStream out(&file);
	out << QString("bcount : %1\n").arg(bcount);
	out << QString("increase texsize : %1 %2 %3\n").	\
	             arg(dblk).arg(rblk).arg(cblk);
      }

      QMessageBox::critical(0,
			    QString("bcount : %1").arg(bcount),
			    QString("increase texsize : %1 %2 %3").	\
			    arg(dblk).arg(rblk).arg(cblk));
      QMessageBox::critical(0,
			    QString("bcount : %1").arg(bcount),
			    QString("bcount : %1").arg(bcount));

    }


  m_compactdepth= dblk * cbsz;
  m_compactrows = rblk * cbsz;
  m_compactcols = cblk * cbsz;

  //----------
  qint64 csize = m_compactdepth;
  csize *= m_compactrows;
  csize *= m_compactcols;
  qint64 dsize = csize;
  if (m_data) delete [] m_data;
  m_data = new uchar[m_nVolumes*dsize];
  memset(m_data, 0, dsize);
  //----------

  uchar *tmp = new uchar[cbsz*cbsz*cbsz];

  for (int l=0; l<=m_minLevel; l++)
    m_cblkno[l] = -1;

  for(int b=0; b<blocklist.count(); b++)
    {
      int ib,jb,kb;
      ib = floor(blocklist[b].x+0.5f);
      jb = floor(blocklist[b].y+0.5f);
      kb = floor(blocklist[b].z+0.5f);
      qint64 bidx = ib*m_rblocks*m_cblocks + jb*m_cblocks + kb;

//      if (bidx >= m_nblocks)
//	QMessageBox::information(0, "", QString("bidx>nblocks : %1 %2 %3").arg(ib).arg(jb).arg(kb));

      if (bidx < m_nblocks)
	loadBlock(slevel,
		  bidx, ib, jb, kb,
		  rblk, cblk,
		  tmp);
    } // b

  delete [] tmp;
}

void
BlockReader::loadAllValidBlocks(int slevel)
{
  if (m_currLod != slevel)
    {
      setLodLevel(slevel);
      int bsz = m_blockSize+m_overlap;
      float frc = (1024.0f*1024.0f)/(m_nVolumes*bsz*bsz*bsz);
      m_maxCacheSize = m_memorySize*frc; // blocks
    }

  memset(m_auxdata, 0, 4*m_nblocks);
  memcpy(m_blockLevel, m_savedblockLevel, m_nblocks);

  int nb = 0;
  for(int b=0; b<m_nblocks; b++)
    if (m_savedblockLevel[b] == 0)
      {
	nb++;
	m_blockLevel[b] = slevel;
	m_auxdata[4*b+3] = 100; // to signify that the block is valid
      }

  m_maxlodLevel = slevel;

  int cbsz = m_blockSize+m_overlap;
  int maxb = m_maxtexSize/cbsz;
  int bcount = nb;
  int cblk = qMin(maxb, bcount);
  int rblk = qMin(maxb, qMax(1, bcount/cblk + (bcount%cblk>0)));
  int dblk = qMin(maxb, qMax(1, bcount/(cblk*rblk) + (bcount%(cblk*rblk) > 0)));

  //QMessageBox::information(0, QString("bcount : %1").arg(bcount), QString("bcount : %1").arg(bcount));

  if (bcount > dblk*rblk*cblk)
    {
      QString flnm = m_fileName[0] + ".error";
      
      {
	QFile file(flnm);
	if (!file.open(QIODevice::WriteOnly |
		     QIODevice::Text))
	  {
	    QMessageBox::critical(0, "Error", "Cannot write to file");
	    return;
	  }

	QTextStream out(&file);
	out << QString("bcount : %1\n").arg(bcount);
	out << QString("increase texsize : %1 %2 %3\n").	\
	             arg(dblk).arg(rblk).arg(cblk);
      }

      QMessageBox::critical(0,
			    QString("Load All Valid Blocks : bcount : %1").arg(bcount),
			    QString("increase texsize : %1 %2 %3").	\
			    arg(dblk).arg(rblk).arg(cblk));
      QMessageBox::critical(0,
			    QString("bcount : %1").arg(bcount),
			    QString("bcount : %1").arg(bcount));

    }


  m_compactdepth= dblk * cbsz;
  m_compactrows = rblk * cbsz;
  m_compactcols = cblk * cbsz;

  //----------
  qint64 csize = m_compactdepth;
  csize *= m_compactrows;
  csize *= m_compactcols;
  qint64 dsize = csize;
  if (m_data) delete [] m_data;
  m_data = new uchar[m_nVolumes*dsize];
  memset(m_data, 0, dsize);
  //----------      

  uchar *tmp = new uchar[cbsz*cbsz*cbsz];

  for (int l=0; l<=m_minLevel; l++)
    m_cblkno[l] = -1;

  for(int bidx=0; bidx<m_nblocks; bidx++)
    {
      if (m_blockLevel[bidx] == slevel)
	{
	  int ib,jb,kb;
	  ib = bidx/(m_rblocks*m_cblocks);
	  jb = (bidx - ib*m_rblocks*m_cblocks)/m_cblocks;
	  kb = bidx - ib*m_rblocks*m_cblocks - jb*m_cblocks;
	  
	  loadBlock(slevel,
		    bidx, ib, jb, kb,
		    rblk, cblk,
		    tmp);
	} // if
    } // bidx

  delete [] tmp;
}

void
BlockReader::loadBlock(int slevel,
		       int bidx, int ib, int jb, int kb,
		       int rblk, int cblk,
		       uchar *tmp)
{
  int cbsz = m_blockSize+m_overlap;
  hsize_t offset[3];

  offset[0] = ib*cbsz;
  offset[1] = jb*cbsz;
  offset[2] = kb*cbsz;

  qint64 mapIdx = slevel*m_nblocks0 + bidx;
  if (m_lru.contains(mapIdx))
    {
      if (m_auxdata[4*bidx+3] >= 200)
	QMessageBox::information(0, "", "Take care before removing this block - we need this block");
      
      m_lru.removeOne(mapIdx);
    }
      
  uchar *block;
  if (m_blockCache.contains(mapIdx))
    memcpy(m_tmpV, m_blockCache[mapIdx], m_nVolumes*cbsz*cbsz*cbsz);
  else
    {
      for(int nf=0; nf<m_nVolumes; nf++)
	{
	  readBlock(nf,
		    slevel,
		    offset,
		    tmp);
	  for(int v=0; v<cbsz*cbsz*cbsz; v++)
	    m_tmpV[m_nVolumes*v+nf] = tmp[v];
	}
      
      if (m_blockCache.size() < m_maxCacheSize)	      
	block = new uchar[m_nVolumes*cbsz*cbsz*cbsz];
      else 
	{
	  block = m_blockCache[m_lru[0]];
	  m_blockCache.remove(m_lru[0]);
	  m_lru.removeFirst();
	}
      
      memcpy(block, m_tmpV, m_nVolumes*cbsz*cbsz*cbsz);
      m_blockCache.insert(mapIdx, block);
    }
  m_lru.append(mapIdx);
  
  
  //- update texture block indices
  m_cblkno[0] ++;
  
  copyBlock(slevel, bidx, rblk, cblk);
}

void
BlockReader::readBlock(int nf,
		       int slevel,
		       hsize_t *offset,
		       uchar *tmp)
{
  int cbsz = m_blockSize+m_overlap;

  hsize_t count[3];
  count[0] = qMin(cbsz, m_odepth- (int)offset[0]-1);
  count[1] = qMin(cbsz, m_orows - (int)offset[1]-1);
  count[2] = qMin(cbsz, m_ocols - (int)offset[2]-1);

  hsize_t mdim[3];
  mdim[0] = mdim[1] = mdim[2] = cbsz;
  hsize_t memoffset[3];
  memoffset[0] = 0;
  memoffset[1] = 0;
  memoffset[2] = 0;
  DataSpace memspace( 3, mdim );
  memspace.selectHyperslab(H5S_SELECT_SET,
			   count,
			   memoffset);

  DataSpace dataspace = m_hdf5Dataset[nf][slevel].getSpace();
  dataspace.selectHyperslab(H5S_SELECT_SET,
			    count,
			    offset);
		    
  m_hdf5Dataset[nf][slevel].read(tmp,
				 m_dataType,
				 memspace,
				 dataspace);
}

void
BlockReader::copyBlock(int slevel, int bidx,
		       int nrows, int ncols)
{
  int v[3];
  v[0] = m_cblkno[0]/(nrows*ncols);
  v[1] =(m_cblkno[0] - v[0]*nrows*ncols)/ncols;
  v[2] = m_cblkno[0] - v[0]*nrows*ncols - v[1]*ncols;
  
  int cbsz = m_blockSize+m_overlap;
  int vistart = v[0]*cbsz;
  int vjstart = v[1]*cbsz;
  int vkstart = v[2]*cbsz;
  
  qint64 csize = m_compactdepth;
  csize *= m_compactrows;
  csize *= m_compactcols;

  int idx = 0;
  for(int i=0; i<cbsz; i++)
    for(int j=0; j<cbsz; j++)
      for(int k=0; k<cbsz; k++)
	{
	  qint64 tidx = ((vistart+i)*nrows*ncols*cbsz*cbsz +
			 (vjstart+j)*ncols*cbsz +
			 (vkstart+k));	  
	  if (tidx >= csize)
	    {
	      int nb0 = 0;
	      for(int b=0; b<m_nblocks; b++)
		if (m_blockLevel[b] == slevel)
		    nb0++;

	      int nb1 = 0;
	      for(int b=0; b<m_nblocks; b++)
		if (m_savedblockLevel[b] == 0)
		    nb1++;

	      QMessageBox::information(0, "error - copy block",
				       QString("%1\n %2 %3\n %4 %5 %6\n %7 %8 %9\n %10 %11").\
				       arg(tidx).arg(nb0).arg(nb1).\
				       arg(m_compactdepth).	   \
				       arg(m_compactrows).	   \
				       arg(m_compactcols).	   \
				       arg(vistart).arg(vjstart).arg(vkstart).\
				       arg(nrows).arg(ncols));
	    }

	  for(int nf=0; nf<m_nVolumes; nf++)
	    m_data[m_nVolumes*tidx+nf] = m_tmpV[m_nVolumes*idx+nf];

	  idx ++;
	}

  m_auxdata[4*bidx+0] = v[2];
  m_auxdata[4*bidx+1] = v[1];
  m_auxdata[4*bidx+2] = v[0];
  m_auxdata[4*bidx+3] = slevel+200;
}

uchar*
BlockReader::wholeVolume(int slevel)
{
  if (m_data)
    delete [] m_data;
  m_data = new uchar[m_depth*m_rows*m_cols];

  m_hdf5Dataset[0][slevel].read(m_data, m_dataType);

  return m_data;
}

uchar*
BlockReader::loadSingleBlock(int ib, int jb, int kb,
			     int slevel)
{
  int bidx = ib*m_rblocks*m_cblocks + jb*m_cblocks + kb;
  if (m_savedblockLevel[bidx] == 255)
    return 0;

  if (m_currLod != slevel)
    {
      setLodLevel(slevel);
      int bsz = m_blockSize+m_overlap;
      float frc = (1024.0f*1024.0f)/(m_nVolumes*bsz*bsz*bsz);
      m_maxCacheSize = m_memorySize*frc; // blocks
    }

  m_maxlodLevel = slevel;
  
  if (bidx >= m_nblocks)
    {
      QMessageBox::information(0, "", QString("%1 %2 %3").arg(ib).arg(jb).arg(kb));
    }

  int cbsz = m_blockSize+m_overlap;
  hsize_t offset[3];
  hsize_t dimsm[3], count[3], memoffset[3];
  dimsm[0] = dimsm[1] = dimsm[2] = cbsz;
  count[0] = count[1] = count[2] = cbsz;
  memoffset[0] = memoffset[1] = memoffset[2] = 0;
  DataSpace memspace( 3, dimsm );

  offset[0] = ib*cbsz;
  offset[1] = jb*cbsz;
  offset[2] = kb*cbsz;

  qint64 mapIdx = slevel*m_nblocks0 + bidx;
  if (m_lru.contains(mapIdx))
    m_lru.removeOne(mapIdx);
      
  if (! m_blockCache.contains(mapIdx))
    {
      uchar *tmp = new uchar[cbsz*cbsz*cbsz];
      for(int nf=0; nf<m_nVolumes; nf++)
	{
	  readBlock(nf,
		    slevel,
		    offset,
		    tmp);
	  for(int v=0; v<cbsz*cbsz*cbsz; v++)
	    m_tmpV[m_nVolumes*v+nf] = tmp[v];
	}
      
      uchar *block;
      if (m_blockCache.size() < m_maxCacheSize)	      
	block = new uchar[cbsz*cbsz*cbsz];
      else 
	{
	  block = m_blockCache[m_lru[0]];
	  m_blockCache.remove(m_lru[0]);
	  m_lru.removeFirst();
	}
      
      memcpy(block, m_tmpV, m_nVolumes*cbsz*cbsz*cbsz);
      m_blockCache.insert(mapIdx, block);
      delete [] tmp;
    }
  m_lru.append(mapIdx);

  return m_blockCache[mapIdx];
}

void
BlockReader::generateHistograms(int nf, int lod)
{
  m_1dHistogram[nf] = new int[256];
  memset(m_1dHistogram[nf], 0, 256*4);

  m_2dHistogram[nf] = new int[256*256];
  memset(m_2dHistogram[nf], 0, 256*256*4);

  float *flhist1D = new float[256];
  memset(flhist1D, 0, 256*4);
  float *flhist2D = new float[256*256];
  memset(flhist2D, 0, 256*256*4);

  setLodLevel(lod);

  int cbsz = m_blockSize+m_overlap;
  uchar *tmp = new uchar[cbsz*cbsz*cbsz];

  QProgressDialog progressBar(QString("Generating histogram from LOD(%1)").arg(lod),
			      "", 0, m_nblocks, 0);
  progressBar.setCancelButton(0);

  int o2 = m_overlap/2;

  int bn = 0;
  for(int d=0; d<m_dblocks; d++)
    for(int r=0; r<m_rblocks; r++)
      for(int c=0; c<m_cblocks; c++)
	{
	  progressBar.setValue(bn);
	  bn++;

	  int dmin = o2;
	  int rmin = o2;
	  int cmin = o2;
	  int dmax = qMin(cbsz-o2, m_odepth-d*cbsz-o2);
	  int rmax = qMin(cbsz-o2, m_orows -r*cbsz-o2);
	  int cmax = qMin(cbsz-o2, m_ocols -c*cbsz-o2);

	  hsize_t offset[3];
	  offset[0] = d*cbsz;
	  offset[1] = r*cbsz;
	  offset[2] = c*cbsz;
	  readBlock(nf,
		    lod,
		    offset,
		    tmp);

	  for(int dm=dmin; dm<dmax; dm++)
	    for(int rm=rmin; rm<rmax; rm++)
	      for(int cm=cmin; cm<cmax; cm++)
		{
		  int idx = dm*cbsz*cbsz + rm*cbsz + cm;
		  flhist1D[tmp[idx]]++;		  
		}

	  for(int dm=dmin; dm<dmax; dm++)
	    for(int rm=rmin; rm<rmax; rm++)
	      for(int cm=cmin; cm<cmax; cm++)
		{
		  int gx = tmp[(dm+1)*cbsz*cbsz + rm*cbsz + cm] -
		           tmp[(dm-1)*cbsz*cbsz + rm*cbsz + cm];
		  int gy = tmp[dm*cbsz*cbsz + (rm+1)*cbsz + cm] -
		           tmp[dm*cbsz*cbsz + (rm-1)*cbsz + cm];
		  int gz = tmp[dm*cbsz*cbsz + rm*cbsz + (cm+1)] -
		           tmp[dm*cbsz*cbsz + rm*cbsz + (cm-1)];

		  int g = sqrtf(gx*gx+gy*gy+gz*gz);
//		  int g = qAbs(gy);
		  g = qBound(0, g, 255);

		  int v = tmp[dm*cbsz*cbsz + rm*cbsz + cm];

		  flhist2D[g*256 + v]++;
		}
	}

  delete [] tmp;

  generateHistograms(flhist1D, flhist2D,
		     m_1dHistogram[nf], m_2dHistogram[nf]);
  delete [] flhist1D;
  delete [] flhist2D;

  progressBar.setValue(m_nblocks);
}

void
BlockReader::generateHistograms(float *flhist1D,
				float *flhist2D,
				int* hist1D,
				int* hist2D)
{
  int i;

  // generate 1d histogram
  float maxf, minf, mlen;
  maxf = -1;
  for (i=1; i<256; i++)
    maxf = qMax(maxf,flhist1D[i]);

  if (maxf > 0)
    {
      for (i=0; i<256; i++)
	hist1D[i] = (int)(255*flhist1D[i]/maxf);
      hist1D[0] = qMin(hist1D[0],255);
    }


  // generate 2d histogram
  for (i=0; i<256*256; i++)
    if (flhist2D[i] > 1)
      flhist2D[i] = log(flhist2D[i]);
    
  maxf = -1.0f;
  minf = 1000000.0f;
  for (i=0; i<256*256; i++)
    {
      maxf = qMax(maxf,flhist2D[i]);
      minf = qMin(minf,flhist2D[i]);
    }

  mlen = maxf-minf;
  if (mlen > 0)
    {
      for (i=0; i<256*256; i++)
	hist2D[i] = (int)(255*(flhist2D[i]-minf)/mlen);
    }
}

Q_EXPORT_PLUGIN2(blockreader, BlockReader);

